module functions_mod
USE parameters_mod
USE parameters_energy_mod
IMPLICIT NONE
CONTAINS

FUNCTION spayer(skinkN,skink,sspay,aime)
INTEGER, intent(in)   ::skinkN
DOUBLE PRECISION,    intent(in)   ::skink(skinkN),sspay(skinkN),aime
DOUBLE PRECISION                  :: spayer,spayer_temp(1),aime_pass(1)

aime_pass(1)=aime
CALL spline_linear_val2(skinkN,skink,sspay,aime_pass(1),spayer_temp,1)
spayer=max(spayer_temp(1),.02D0)

end FUNCTION spayer



FUNCTION c_upper_fnc2(labor,effective_wage,capital,kprime,kdim,tr,r)
INTEGER, intent(in)   ::kdim
DOUBLE PRECISION,    intent(in)   ::kprime(kdim),capital,labor,effective_wage
DOUBLE PRECISION,    intent(in)   ::tr,r
DOUBLE PRECISION                  :: c_upper_fnc2(kdim), c_upper_temp(kdim)

INTEGER               :: i


do i=1,kdim

    c_upper_temp(i)=workcon(labor,capital,kprime(i),tr,r,effective_wage)
end do
c_upper_fnc2=c_upper_temp
END FUNCTION c_upper_fnc2

FUNCTION c_upper_fnc(prod,capital,kprime,kdim,tr,w,r)
IMPLICIT NONE
INTEGER, intent(in)   ::kdim
DOUBLE PRECISION,    intent(in)   ::kprime(kdim),capital
DOUBLE PRECISION,    intent(in)   ::tr,w,r,prod
DOUBLE PRECISION                  :: c_upper_fnc(kdim), c_upper_temp(kdim)

INTEGER               :: i

do i=1,kdim
    c_upper_temp(i)=workcon(0.90D0,capital,kprime(i),tr,r,w*prod)
end do
c_upper_fnc=c_upper_temp
END FUNCTION c_upper_fnc




FUNCTION consumption_tilda(con,eng,adl_equiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: con,eng,adl_equiv
DOUBLE PRECISION                  :: consumption_tilda
consumption_tilda=(con**gamma1*(eng-ebar)**(1.0D0-gamma1))/adl_equiv
END FUNCTION consumption_tilda



FUNCTION workcon(labor,capital,kprime,tr,r,w)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital,kprime,tr,r,labor,w
DOUBLE PRECISION                  :: workcon,ss_taxable

ss_taxable=min(labor*w,ss_cap)
workcon=(1+r)*(capital+tr)+labor*w-ss_taxable*tauss-kprime-all_taxer(r*(tr+capital),w*labor)
END FUNCTION workcon


FUNCTION all_taxer(capital_income,labor_income)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income
DOUBLE PRECISION                  :: all_taxer,ypp
ypp=labor_income-.5D0*tauss*min(ss_cap,labor_income)
all_taxer=lab_taxer(ypp)+capital_income*lambdak
end FUNCTION all_taxer

FUNCTION all_taxer2(capital_income,labor_income,lambda0_test)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital_income,labor_income,lambda0_test
DOUBLE PRECISION                  :: all_taxer2,ypp
ypp=labor_income-.5D0*tauss*min(ss_cap,labor_income)
if(ypp==0.0D0) then
    all_taxer2=0.0D0+capital_income*lambdak
else
all_taxer2=(1.0D0-lambda0_test*((ypp/avg_earnings_scale)**(-lambda1)))*ypp+capital_income*lambdak
end if
end FUNCTION all_taxer2


FUNCTION lab_taxer(income)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: income
DOUBLE PRECISION                  :: lab_taxer
if(income==0.0D0) then
    lab_taxer=0.0D0
else
lab_taxer=(1.0D0-lambda0*((income/avg_earnings_scale)**(-lambda1)))*income
end if

end FUNCTION lab_taxer


FUNCTION ret_taxer(capital,kprime,tr,r,ss_pass,avg_earnings)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital,kprime,tr,r,ss_pass,avg_earnings
DOUBLE PRECISION                  :: ret_taxer,tax_ss_income
    tax_ss_income=.0D0*ss_pass

    ret_taxer=lab_taxer(tax_ss_income)

END FUNCTION ret_taxer

FUNCTION retcon(capital,kprime,tr,r,ss_pass,avg_earnings)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: capital,kprime,tr,r,ss_pass,avg_earnings
DOUBLE PRECISION                  :: retcon,tax_ss_income
    tax_ss_income=.0D0*ss_pass

    retcon=(1+r)*(capital+tr)+ss_pass-lambdak*(r*(tr+capital))-kprime-lab_taxer(tax_ss_income)

END FUNCTION retcon



FUNCTION retcon_ar(dim,capital,kprime,tr,r,ss_pass,avg_earnings)
IMPLICIT NONE
INTEGER,    intent(in)   :: dim
DOUBLE PRECISION,    intent(in)   :: capital,kprime(dim),tr,r,ss_pass,avg_earnings
DOUBLE PRECISION                  :: retcon_ar(dim),tax_ss_income
INTEGER               :: it
tax_ss_income=.0D0*ss_pass

do it=1,dim
    retcon_ar(it)=(1+r)*(capital+tr)+ss_pass-lambdak*(r*(tr+capital))-kprime(it)-lab_taxer(tax_ss_income)
end do
END FUNCTION retcon_ar




FUNCTION utility_last(con,eng,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: con,eng,ad_eqiv
DOUBLE PRECISION                  :: utils_tmp, utility_last,lc
if ( (eng+con)>0.015 .and. eng>ebar) then
    lc=consumption_tilda(con,eng,ad_eqiv)
    utils_tmp =lc**(1.0D0-sigma1)/(1.0D0-sigma1)*ad_eqiv
else
    utils_tmp = -1000000000000000000.0D0
end if

utility_last = utils_tmp
END FUNCTION utility_last


FUNCTION valuer(c,e,labor,discount,vprime,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: c,e,labor,ad_eqiv
DOUBLE PRECISION,    intent(in)   :: discount,vprime
DOUBLE PRECISION                  :: utils_tmp,valuer,lc
    if ( (c+e)>0.015 .and. e>ebar) then
        lc=consumption_tilda(c,e,ad_eqiv)
        if ( labor>0.01D0) then
            utils_tmp = ad_eqiv*(lc**(1.0D0-sigma1)/(1.0D0-sigma1)-chi*labor**(1.0D0+1.0D0/sigma2)/(1.0D0+1.0D0/sigma2))+discount*beta*vprime
        else
            utils_tmp =ad_eqiv*lc**(1.0D0-sigma1)/(1.0D0-sigma1)+discount*beta*vprime
        end if
    else
        utils_tmp = -1000000000000000000.0D0
    end if


    valuer = utils_tmp
END FUNCTION valuer

FUNCTION valuer_ret(c,e,discount,vprime,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: c,e,ad_eqiv
DOUBLE PRECISION,    intent(in)   :: discount
DOUBLE PRECISION                  :: utils_tmp, valuer_ret,vprime,lc
      if ( (c+e)>0.015 .and. e>ebar) then
        lc=consumption_tilda(c,e,ad_eqiv)
        utils_tmp = ad_eqiv*lc**(1.0D0-sigma1)/(1.0D0-sigma1)+beta*discount*vprime
      else
         utils_tmp = -1000000000000000000.0D0
      end if

valuer_ret = utils_tmp
END FUNCTION valuer_ret

FUNCTION valuer_ret1(c,e,discount,vprime,ad_eqiv)
IMPLICIT NONE
DOUBLE PRECISION,    intent(in)   :: c,e,ad_eqiv
DOUBLE PRECISION,    intent(in)   :: discount
DOUBLE PRECISION                  :: utils_tmp, valuer_ret1,vprime,lc
      if ( (c+e)>0.015 .and. e>ebar) then
        lc=consumption_tilda(c,e,ad_eqiv)
        utils_tmp = ad_eqiv*lc**(1.0D0-sigma1)/(1.0D0-sigma1)+beta*discount*vprime
print*,"valuer_ret1_out",lc,c,e,beta,discount,vprime
      else
         utils_tmp = -1000000000000000000.0D0
      end if

valuer_ret1 = utils_tmp
END FUNCTION valuer_ret1

SUBROUTINE energy_root_finder(inputs,misses,choicesN)
!
IMPLICIT NONE
INTEGER, intent(in)   ::choicesN
DOUBLE PRECISION,    intent(out)  :: misses(choicesN)
DOUBLE PRECISION,    intent(in)   :: inputs(choicesN)
DOUBLE PRECISION   :: Ep_guess,E_guess,K1_guess,K2_guess,N1_guess,N2_guess,K_in,N_in,Ec_in,Ae_in,Atfp_in,pe_in,taue_in
K_in=passed_parameters(1)
N_in=passed_parameters(2)
Ec_in=passed_parameters(3)
Atfp_in=passed_parameters(4)
Ae_in=passed_parameters(5)
taue_in=passed_parameters(6)

Ep_guess=inputs(1)
K1_guess=inputs(2)
N1_guess=inputs(3)

K2_guess=K_in-K1_guess
N2_guess=N_in-N1_guess
E_guess=Ae_in*K2_guess**alpha2*N2_guess**(1.0D0-alpha2)

pe_in=theta*Atfp_in*(K1_guess/Ep_guess)**alpha1*(N1_guess/Ep_guess)**(1.0D0-alpha1-theta)-taue_in
if(N1_guess<0.0D0 .or. N2_guess<0.0D0 .or. K1_guess<0.0D0 .or. K2_guess<0.0D0 .or. Ep_guess<0.0D0) then
misses(1)=100.0D0*(max(abs(N1_guess),0.0D0)+max(abs(N2_guess),0.0D0)+max(abs(K1_guess),0.0D0)+max(abs(K2_guess),0.0D0)+max(abs(E_guess),0.0D0)+max(abs(Ep_guess),0.0D0))
misses(2)=misses(1)*10.0D0
misses(3)=misses(1)*5.0D0
else
misses(1)=E_guess-Ep_guess-Ec_in
misses(2)=Atfp_in*alpha1*(N1_guess/K1_guess)**(1.0D0-alpha1-theta)*(Ep_guess/K1_guess)**theta-pe_in*alpha2*Ae_in*(N2_guess/K2_guess)**(1.0D0-alpha2)
misses(3)=Atfp_in*(1.0D0-alpha1-theta)*(K1_guess/N1_guess)**alpha1*(Ep_guess/N1_guess)**theta-pe_in*Ae_in*(1.0D0-alpha2)*(K2_guess/N2_guess)**alpha2
end if

END SUBROUTINE energy_root_finder



FUNCTION energy_share_equation(c_con,totconsumption,age)
IMPLICIT NONE
INTEGER, intent(in)   ::age
DOUBLE PRECISION,    intent(in)   :: c_con,totconsumption
DOUBLE PRECISION   :: e_con,energy_share_equation

e_con=(totconsumption-c_con)/pe

energy_share_equation =&
-consumption_tilda(c_con,e_con,adult_equivalence(age))
END FUNCTION energy_share_equation




DOUBLE PRECISION FUNCTION energy_sharer(ax,bx,cx,func,tol,xmin,totconsumption,age)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    INTEGER, intent(in)   ::age
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,totconsumption
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,totconsumption,age)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        energy_sharer=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,totconsumption,age)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('energy_sharer: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION energy_sharer



DOUBLE PRECISION FUNCTION energy_leveler(ax,bx,cx,func,tol,xmin,K,N)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,K,N
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,K,N)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        energy_leveler=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,K,N)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('energy_leveler: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION energy_leveler


subroutine spline_linear_val2( ndata, tdata, ydata, tval, yval, ndata2 )
      INTEGER :: left, right,i
      INTEGER, intent(in)  ::ndata,ndata2
      DOUBLE PRECISION,  intent(in) :: tdata(ndata),ydata(ndata),tval(ndata2)
      DOUBLE PRECISION,  intent(out) :: yval(ndata2)
      DOUBLE PRECISION  ::  ypval
ypval=0.0D0
yval=0.0D0
do i=1,ndata2
    call rvec_bracket2 ( ndata, tdata, tval(i), left, right )
  ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) )
  yval(i) = ydata(left) +  ( tval(i) - tdata(left) ) * ypval
end do
  return
end SUBROUTINE spline_linear_val2

subroutine rvec_bracket2 ( n, x, xval, left, right )
  implicit none
      INTEGER, intent(in)  ::n
      INTEGER, intent(out)  ::left,right
      DOUBLE PRECISION, intent(in)   ::xval,x(n)
  integer i
  do i = 2, n - 1

    if ( xval < x(i) ) then
      left = i - 1
      right = i
      return
    end if
   end do
  left = n - 1
  right = n
  return
end subroutine rvec_bracket2

END module functions_mod

